Single electron transistor strongly coupled to vibrations: Counting Statistics and 

Fluctuation Theorem 
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Using a simple quantum master equation approach, we calculate the Full Counting Statistics of a 
single electron transistor strongly coupled to vibrations. The Full Counting Statistics contains both 
the statistics of integrated particle and energy currents associated to the transferred electrons and 
phonons. A universal as well as an effective fluctuation theorem are derived for the general case 
where the various reservoir temperatures and chemical potentials are different. The first relates to 
, the entropy production generated in the junction while the second reveals internal information of 

the system. The model recovers Franck-Condon blockade and potential applications to non-invasive 
f~* > ; ' molecular spectroscopy are discussed. 



PACS numbers: 05.60.Gg, 73.23.Hk, 05.70.Ln, 74.25.fg 



The Full Counting Statistics (FCS) of energy and matter exchanges provides a wealth of information about the 
dynamics of multi-terminal nanostructures. A great achievement of the last decade has been to conveniently identify 
universal features in the FCS [1-6]. Roughly speaking, these arc related to the fact that fluctuations of entropy 
production, AjS 1 , satisfy a universal fluctuation theorem (FT) P+Ais/P-AiS = e 1 ■ At steady state, the entropy 
'-^ ■ production in the nanostructure must be balanced by the entropy flow through its terminals [7] , which can be accessed 
[ through the FCS. This implies that entropy production is a measurable quantity and can be expressed as a sum of 

■ the various thermodynamic affinities acting on the system times their associated fluxes [8] . Fostered by the increased 
experimental abilities enabling the counting of single electron transfers [9] and thereby the experimental verification 
of the FT [10], recent works have identified non- universal FT-symmetries when an incomplete monitoring of all the 
terminals is performed. The resulting effective affinities may be used to probe system-specific features [11-15]. 

Nanoscale devices displaying coupling between electronic and vibrational transport have been significantly studied 
in the past, in part due to their importance for thermo-electricity [16-18]. Most studies rely on the weak coupling 
assumption and those who don't use sophisticated self-consistent nonequilibrium Greens functions procedures which 
often restrict the study of the FCS to the first few moments [19-22]. 

In this paper, we calculate the FCS for a single electron transistor (SET) strongly coupled to a phonon bath using 

■ a new particularly simple quantum master equation (QME) approach. We explicitly derive the universal FT in the 
general case when different electronic and phononic temperatures are considered. Effective FTs are also identified and 
both the electron and phonon average currents and Fano factors are analysed. Going beyond heat exchange [23, 24], 
our model provides an efficient probe for non-destructive molecular spectroscopy: While it is known that one can 

C^l , determine phonon energies from electronic transport characteristics only, the simplicity of our model allows to identify 

o ■ parameter regimes where this is possible without inducillg furthcr heating of phonons ' 
m 



I. MODEL 

We consider the Hamiltonian 

H = eSd + e ka cl a c ka + ^ io q a\a q 

ka q 

+ [*fa»dcL +h - c ] +^d^2[h q a q + h.c.} (1) 



describing a SET with two electronic leads (treated perturbatively) and additionally coupled to one or many phononic 
modes (treated non-perturbatively) . Fermionic operators d (cka) annihilate electrons on the dot (lead a <G {L,R}) 
with energies e ( and a q are the bosonic annihilation operators for a phonon with energy u q . The parameters t ka 
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and h q describe electronic tunneling and phonemic absorption amplitudes, respectively. For only a few phonon modes 
q G {1, . . . , Q}, the model may e.g. describe electronic transport through a molecule where the electronic occupation 
couples to molecular vibrations. In particular for a single phonon mode (Q = 1), the model represents a special 
case of the Anderson-Holstein model, which has previously been treated for example in the linear response and weak 
electron-phonon-coupling regimes [25] or with focus on the electronic charge transfer statistics [26, 27]. For many 
different phonon modes or even a continuum, the model may describe more complex molecules or the interaction with 
bulk phonons, respectively. 

Establishing the complete electronic and phononic FCS requires to monitor not only the charge transfer but also 
the emitted and absorbed bosons. In principle, the statistics of the latter can be retrieved by modelling the phonons 
as part of the system, see e.g. Refs. [28-31]. This allows to explore the strong electron-phonon coupling limit, but also 
leads to an infinitely large Hilbert space, which renders the study of the full dynamics tedious and hard to interpret. 
For many different phonon modes or even a continuum the required computational resources make it completely 
infeasible to follow this approach. Here, we therefore aim at an efficient representation with only the dot occupation 
treated as a dynamical variable whilst all electronic and phonon terminals are held at thermal equilibrium states. The 
challenge is to retain the complete FCS from such a reduced model. We would like to emphasize here that even in 
the single phonon mode case the postulated stationarity of the phonons does not imply a weak coupling assumption 
between electrons and phonons. For example, for strongly coupled phonons there might exist an even faster relaxation 
process for the phonons immediately restoring thermal equilibrium. 

After a polaron (Lang-Firsov) transformation [32, 33] H' = e +s He~ s with S = $ d^2 q {-^-a q — ^-a q ), where 



ka q 

V ( h " q n"<' hq n 

ka 



tkadc J ka e ' v 5 / + h.c 



(2) 



the dot energy is renormalized e = e — y\ J , and the electronic tunneling is dressed by exponentials of bosonic an- 
nihilation and creation operators. Expanding these terms demonstrates that each electronic tunneling event may now 
be accompanied by multiple phonon emissions and absorptions. Most important however, we note that a perturbative 
treatment in the electronic tunneling amplitudes t)~ a - valid for small electronic tunneling rates in comparison to the 
electronic reservoir temperatures - still allows for a non-perturbative treatment of the electron-phonon interaction 
(parametrised by h q ). 

We now apply standard techniques (e.g. [2, 34]) to derive a QME containing the FCS of both the emitted/absorbed 
phonons and the electrons having traversed the system. Note that since we are perturbative in the electronic tunneling 
amplitudes, we neglect coherences and also the Kondo effect remains out of reach. For a system-bath decomposition 
of the form H§b = J^i^iBi with system and bath operators Ai and Bi, respectively, this QME requires to calculate 
the bath correlation function CV,(t) = (e +lffBT i?ie -liiBr -B :) ), where the expectation value is taken with respect to the 

stationary reservoir state [35]. In our case, the latter consists of a tensor product pe = Pb^ ® Pb^ ® Pb^ °f different 
equilibrium states characterized by temperatures, and for the electronic leads also by chemical potentials. We choose 
the system coupling operators as Ai tCt = d and Ai^ a = dfi and correspondingly we see that the contribution from 

the phonon bath enters multiplicatively in the bath coupling operators B\ a = ^ fe tk a c\ a exp[— XXz!r a <? — ZT a q)] an< ^ 

B2, a = B\ a . In contrast, the contributions from the two electronic leads enter additively as usual. This leads (see A 

for more details) to a non-standard product form of the correlation function Ce(r) = C"' ol (T)C p h(T), where 
I € {(12), (21)}. For the electronic contribution to the correlation function, we have the usual Fourier decomposition 
C"' el (r) = J j"' cl (uj)e~ iuJT duj with the standard electronic Fourier transforms 

7l2 el ( w ) = r Q (-w)/ Q (-w) , 721 C V) = r «(+ w ) [1 - Za(+W)] , (3) 

where f a (uj) = [ e Pn(u-fj. a ) _|_ 1 d cno tcs the Fermi function of lead a held at inverse temperature /3 a and chemical 
potential [i a . The electronic tunneling rates in Eq. (3) are defined by T a (ui) = 27r^fc K/ca| 2 <5(w — tka)- In the following, 
we will parameterize them by a Lorentzian shape T a {ui) = T a 8 2 / {ui 2 + 5 2 ) with S — > oo characterizing the wide-band 
limit. Finite S may effectively model non-Markovian effects induced by the electronic environment [36]. Combining 
the phonon contributions to the bath correlation function using the Baker-Campbell-Hausdorff formula, we obtain 

C ph (r) = exp M [e~^ (1 + n|) + e+^n* - (1 + 2n|)] ) (4) 
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with n B = ^e^p 11 " 5 — l] denoting the Bose distribution with inverse temperature /3 p h- Wc note that the Fourier 
transform of the combined correlation function 7^(0;) = J Ce(T)e +luJT d,T constitutes the transition rates in the rate 
equation when evaluated at the renormalizcd dot energies ±e. 



II. FLUCTUATION THEOREMS 



A. Entropy Production 

For a finite number of phonon modes, the FT for entropy production may be expressed in terms of the FCS, which 
is detailed in B. This is most accessible for a single phonon mode (Q = 1) at frequency u>i = f2 (we also abbreviate 
rig — >■ 7Jb and hi —> h), where the sum in the exponential of Eq. (4) collapses. Then, one can identify the Fourier 
transform 7™(w) of the correlation function via a simple integral transformation. In particular, it is straightforward 
to see (see A) that it can be decomposed 7"(w) = Xn^oo 7"n( w ) m *° processes associated with the net absorption 
(emission) of n > (n < 0) quanta by the phonon modes 



a 1 \ Q.ol/ n _Jillf.n+2nB) / 1 + "B 



"B 



i/2 

X 



M 2 



xJn ^2^- v /n B (l + n B )J , (5) 

where J n {x) denotes the modified Bessel function of the first kind. First, in the zero-coupling limit |/i| 2 — > 0, 
we recover the correlation functions of the SET (3) as all contributions with n ^ vanish. Second, in the 
zero-phonon-temperature limit tib — > 0, only contributions for absorption by the phonon bath remain with 
lim [l+ns.ln/2 J wWl + n B )) "=° 0-^-) nl j. Third, in the wide-band (5 00 such that T a ko) -> T a ) and 

n B — >0 ,lB Si n. 

infinite bias limit (/i Q — > {—00, 00} so that /a(w) — > {0, 1}), we recover the standard infinite-bias results of the SET. 

Both the Fermi-Dirac distribution and the Bose-Einstein distribution obey the separate relations 
/ a (w) = e-Pai^'-t 1 '') [1 _ f a (u>)] and tib = e~@ phn (1 + n B )- These imply that the correlation functions (5) obey a 
relation of Kubo-Martin-Schwinger (KMS) type involving both electronic and phononic temperatures, 

7i\ + n(-^) = e-^("-^ + " 'e^" fi 72 Q l! _„(+ W ), (6) 

which in the case of equal temperatures /3 p h = Pl = Pr reduces to the conventional KMS condition [37]. Evaluated 
at the renormalized dot level e = e — \h\ 2 /il, these enter the Liouvillian C of the QME describing the dot dynamics: 
p = Cp, where the part acting on the populations only reads 



a£{L,R} 5 



- — 00 



Here, jf 2 n = 7?2 n (~^) denotes the transition rate from an empty to a filled dot due to an electron jumping in from 
lead a whilst simultaneously triggering the absorption of n a quanta by the phonons. Similarly, 7^ _„ a = J21 -n (+^) 
denotes the rate of the inverse process. Effectively, this corresponds to transport through an electronic system with 
infinitely many reservoirs with shifted chemical potentials, and the energetics of heat transfers is easily accessible on 
the level of single trajectories, see Fig. 1. In order to monitor the full entropy flow through the junctions, it is in 
general necessary to monitor both energy and particle flows. The electronic counting fields Xa count the net number 
of electrons n^f jumping out of the SET to the lead a, while the phonon counting fields £ Q count the corresponding 
net energy transfercd to the phonons ~ n Q Sl. Together they can be used to determine the energy-particle FCS of 
our model, sec B, i.e., to calculate the time-dependent probability distribution P («> (l) (R)(t) of electronic matter 

and associated phonon energy transfers. 

Eq. (7) is a key result of our paper which clearly shows that when phonons are kept thermally equilibrated, a compact 
description of the dynamics can be established which keeps the full information about the electron and phonon counting 
statistics also in the strong electron-phonon coupling regime. In addition to thcrmalized phonons, the perturbative 
treatment of the electronic tunneling requires that /3 Q r Q <C 1. Our approach is consistent thermodynamically since 
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FIG. 1: Sketch of the heat transfer energetics at the trajectory level for the single-phonon-mode (Q = 1) case and for two 
sample elementary processes (black solid with specific heat flows into the system noted). Filled regions left and right denote 
Fermi functions of left and right electronic leads whereas the sketched oscillator denotes the phonon reservoir, respectively. 
During the electronic tunneling events (wavy lines), energy from the phonons may assist electrons far from the dot level e to 
participate in transport. In the left trajectory, an electron jumps from the left lead to the initially empty dot via the emission 
of two quanta by the phonons. In the right trajectory, it leaves the dot to the right lead via the emission of two quanta by 
the phonons. For trajectories with identical initial and final SET state (the change of SET entropy is zero), the total entropy 
production is minus the entropy flow and is given by AiS = -Q { e f> /T L - /T R - (Q p ^ + Q$)/T ph , with transferred 
electronic heat = — (enffi — e ph ) + A t a n if' > an d associated phonon heat Q p °^ — —e p °^. For example, the combination of 
both trajectories shown yields with n^f = — 1, n *^ = +1, = — 2fi, and = — 3f2 altogether a total entropy production 
of A;S* = -(£ - 2J2 - hl)/T l - 5fi/T ph + (i + SQ - fi R )/T R . 



as shown in B, the following universal FT can be derived (also for multiple phonon modes) 

+"oi .+« P h .+e ph \ t^oo , (R) , (L) . (R) (H) 
P-tf>,-$>,-4»®) ~ 61 +AL ^ +ARe ^ ' (8) 

where the affinities corresponding to the phonon energy and electronic fluxes are respectively given by 

A a = /3 ph - p a , A = (fa- p L )e + (l3 L fi L - PrUr). (9) 

We now turn to the discussion of Eq. (8). Each elementary transfer process described by QME (7) involves an 
electron transfer between the SET and a lead a coupled to an energy transfer with the phonons. Since energy is 
conserved during such transfers, the energy that when rv^' > leaves the SET, en^\ is equal to the energy sent 
to the phonons, , plus the energy sent to lead a, eh^j — . As a result the heat entering the SET from the 
phonons is given by = — e^" and the heat entering the SET from lead a by Q^' = — (en^' — ) + ^ a n^ ■ 
For a single phonon mode, = n a Cl, while for multiple modes, = '^2 q bJ q n a,qi where n a and n a , q denote the 
number of quanta with frequency f2 and uj q that are emitted from the phonons (see B for details). The entropy flow 
associated to this process (corresponding to minus the change in the entropy of the phonon and that of the lead a) 
will be given by A e S^ = /3 a Q^' + PphQ p ^ [2]. At steady state, in average and in a large deviation sense (and in 
a strict sense for trajectories connecting identical initial and final SET states), entropy production is equal to minus 
the entropy flow, i.e. A^S = — J2 a {PaQei + /^phQph )■ Also, the total number of transferred electrons has to be 

conserved + =0. As a result, we easily verify that the entropy production A{S becomes equal to Eq. (8). 

The FT (8) thus states that trajectories with a positive entropy production occur with larger probabilities than 
the inverse trajectories, i.e., the expectation value of the entropy production is always positive as predicted by 
the second law. This holds far from equilibrium and non-perturbatively in the electron-phonon coupling strength. 
Single trajectories - occurring with an exponentially suppressed probability - however may have a negative entropy 
production. For example, when fj,£ = /j,r, /3l = Pr = Peh an d Pph > Pel, the combination of the two elementary 
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processes depicted in Fig. 1 with n^' = +1, = — 2fi, and = — 3S1 would yield a trajectory with negative 
production. 

As further consistency tests of Eq. (8) we mention that when (3 C \ = /3r = /?£, we get A = /3 e \(jiL — Hr) and 
= Pph — Pel- When furthermore ft = /3 e i = /5 p h, = and the entropy production becomes identical to that 

of an isolated junction, i.e., Ai5 — > /3Vbi as ?i<ii^ with bias voltage Vbias = Mi — Mft- All affinities obviously vanish at 
equilibrium where (3l = = /? p h and = /I/?. 

Finally, we mention that when phonons were treated as part of the system, the derivation of the FT would be 
standard due to the additivity of the electronic leads: It would not involve entropy flow across a strongly-coupled 
terminal and thus contain no phonon-related affinity. 



B. Incomplete Fluctuation Theorem 



Returning to the general case but disregarding now the phonon heat counting, which is technically performed by 
setting £ Q = 0, an incomplete FT [11-15] is obtained for the probability P (R)(t) of electrons having crossed the 
junction after time t 



In 



t-Hx> (R) 



hi 



7l2(-g)72l( + g) 

7f 2 (-?)72 L i(+e) 



(10) 



where the effective affinity a is not universal anymore (unless (3l = /3r = /3 p h) and can be evaluated numerically, see 
Fig. 2. Consequently, for an incomplete FT ■ a cannot be associated anymore with the total entropy production. 




electron vs. phonon temperature TJT 



FIG. 2: Effective affinity a, Eq. (10), for different values of the coupling strength A = |/i| 2 /n 2 in the single-mode case (solid 
lines) or A = Jo in the continuum case (dashed curves). When T e \ = T p h, the FT is independent of the coupling strength 
(intersection point). For vanishing coupling strength, we recover a = /3 e i(/xz — /j.r) (dotted thin curve). Other Parameters: 
Tl = r fl , /3 p h/ii = +5, PphfJ-B. = —5, S = ion, /3 p h£ = 0, and /3 p h£l = I3 ph iu c = 1. 

However, the non-universal a does now provide additional information about the system: In the example shown in 
Fig. 2, it may e. g. be used as an effective thermometer for the electron temperature T e \ at fixed phonon temperature 
Tph- Note that fiat (S — > oo) electronic tunneling rates cancel in Eq. (10). 
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For a quantum dot coupled to bulk phonons (continuum of modes), obtaining the resolved FCS of different modes 
is realistically impossible in an experiment. Therefore, monitoring only the electronic statistics will also yield an 
incomplete FT. Then, the sum in Eq. (4) can be converted into an integral by introducing the spectral coupling 
density J(cj) = J2 q \h q \ <5(w — ui q ). The common choice of an Ohmic parameterization [33] J{oS) = Jowe~"/ LJc with 
the dimcnsionless coupling strength Jo and cutoff frequency uj c enables one to calculate the phonon contribution to 
the correlation function (4) analytically (we use an overbar to denote the continuum case) 



Cph(r) 



/3ph"c ) \ /3ph"c 



r 2 (^)(i + i™ c ) 



(11) 



where T{x) denotes the Y- function. It is straightforward to verify a separate KMS condition in the time-domain 
Cph( r ) = C p h(—T — i/3 p h), which implies 7 p h(w) = e +/3ph "7 P h(— w) in the frequency domain. The full Fourier trans- 
forms of the correlation functions are therefore given by a convolution integral of the separate Fourier transforms 
7"(f2) = J 7"' ol (w)7 p h(ri — oj)doj. Using the separate electronic and phonon KMS conditions demonstrates that 
for equal temperatures, the conventional KMS condition must also hold in the continuum case. In addition, in the 
wide-band and infinite bias limits - where the electronic Fourier transforms (3) become constants (over a sufficiently 
wide range) - the convolution representation implies that the corresponding electronic currents arc unaffected by the 
additional phonon reservoir. When only electrons are counted, these enter the QME as 



-7i Q 2 (-?) +7 2 a i(+e)e+ i *° 
-7i a 2 (-e>- ix " -7 2 "(+e) 



(12) 



Again it suffices to consider only a single electronic counting field to obtain an incomplete FT as in Eq. (10) with 
7" — > 7", where the renormalized dot level now becomes e = e — Jql>j c . When we choose different temperatures for 
electrons and phonons, the affinity a in the FT may be evaluated numerically, and displays a similar dependence as 
in the single-mode case provided Jq = h 2 /H 2 and lu c = fi, cf. dashed lines in Fig. 2 and the discussion below. 



III. COUNTING STATISTICS 

The cumulants of the FCS such as electronic current I c \ = j| lim (n^\ and noise S e i = 

22 lim ((^n^ ^ — (^cf ^ ^ can be easily calculated using Eqs. (7) and (12) - see e.g. the appendix. We plot 

current and Fano factor F = S/\I\ in Fig. 3. At low electronic temperatures, and in the single-mode case, the phonon 
frequency can be deduced from the steps in the electronic current. Generally, for strong SET-phonon coupling one 
observes a suppression of the finite-bias current in comparison to the uncoupled case - known as Franck-Condon 
blockade [28], which normally requires to treat phonons as dynamical degrees of freedom. Remarkably, our simple 
2x2 rate equations even recover that the giant Fano factors typical for Franck-Condon blockade can only be observed 
when electron and phonon temperatures are chosen very different, i.e., when the reservoirs are unequilibrated. This 
behavior is found in the continuum phonon case (cf. dashed curves in Fig. 3) as well as in the large bias regime 
(/l(cj) 1 and fn(uj) — > but keeping the energy- dependent tunneling rates). In both wide-band and infinite bias 
regime we simply recover the conventional current of the SET (not shown). 

A common problem in molecular transport spectroscopy is that the electronic current used to probe the phonon fre- 
quency will, at high bias, induce vibronic excitations, which may eventually lead to the destruction of the molecule [38]. 
The low-dimensional model proposed here is quite useful to find regimes allowing electronic spectroscopy (finite bias) 
simultaneously with a net negative phonon energy current (implying a lowering of the effective molecular temper- 
ature), see the red solid curve in Fig. 4. There, we evaluate the statistics of total phonon heat emission in the 
single-mode case from Liouvillian (7) by not making any difference between right and left associated electronic jumps, 
see the appendix. The particular parametrization in Fig. 4 demonstrates a stationary configuration with only slightly 
different electronic and phonon temperatures T e \ sa 0.86T p h at finite bias (intersection point of the solid red curve with 
1 = marked by the green vertical line), at which the phonon energy may be probed without inducing further heating. 
More generally, determining the phonon frequency by electronic spectroscopy in a non-destructive manner requires 
electron and phonon temperatures where, within a finite voltage range, the phonon heat current is non-positive. A 
major advantage of our model is that this range can be estimated without having to consider hundreds or thousands 
of phonon states. 
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bias voltage V,. IQ. 



FIG. 3: Electronic current and Fano factor (inset) for the single-mode (solid) and continuum (dashed) cases for different 
temperature configurations (colored curves) and different coupling strengths (brown solid and black dashed). Low electronic 
temperatures are required to resolve the single-mode nature of the phonon bath (steplike solid curves vs. smooth dashed 
curves). The reduction of even the large-bias currents is a consequence of a constant and finite width of the electronic tunneling 
rates. Other Parameters: Yl = T_r = T, 5 = 10f2, e/fi = 0, ii = uj c , and for the bold colored curves A = |/i| 2 /f2 2 = Jo = 5. 



IV. SUMMARY 

We studied the combined Full Counting Statistics of electrons and phonons for a special case of the Anderson- 
Holstein model using an extremely efficient two-dimensional rate equation for the electronic populations which treats 
the phonon mode as a strongly-coupled reservoir. Despite its simplicity, our approach was able to predict strong 
signatures of Franck-Condon blockade in the electronic current and noise. Furthermore, we proved its thermodynamic 
consistency by explicitly deriving a universal fluctuation theorem for the entropy production and identifying the 
three thermodynamic affinities characterizing our setup and imposing universal symmetries on the Full Counting 
Statistics. Our results hold for arbitrary electron-phonon coupling and also for multiple phonon modes as long as 
they are held at thermal equilibrium. In addition, we found a non-universal fluctuation theorem referring only to 
the electronic counting statistics, which may be used to probe the phonon temperature. The universal fluctuation 
theorem is recovered in this case when all terminal temperatures are set equal, which implies that the energy flow 
between the system and the phonons vanishes. Due to its small dimensionality, our model may prove useful to identify 
appropriate regimes for performing non-destructive molecular spectroscopy. Furthermore, we expect our approach to 
be also applicable to electronic transport through more complex structures. 
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FIG. 4: Total energy current into the phonon bath versus electronic temperature T e i/T p h for different phonon temperatures 
/3 p hfi = 0.1 (solid), /3phfi = 1-0 (dashed) and vanishing Vbias^ = (black), intermediate Vbias^ = 15 (red), and high Vbias^ = 30 
(blue) electronic bias voltages. Other Parameters: Yl = T_r, S = 10f2, e/Q = 0, f2/3 p h = 1, \h\ 2 /D. 2 = 5. 
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Appendix A: Bath correlation functions for discrete phonon modes 



Obviously, the phonon correlation function Eq. (4) can be written as a product of single-mode correlation functions 

Q 

Cph(r) = [] Cph(r) , (Al) 

9=1 



where Q denotes the number of different phonon modes. Each factor can be formally expanded in the variables e 



±icj r 



q J mlm'l 



= E C ?n m ', P ^r). (A2) 

m,m'— 

Now, the key observation is that the terms in the sum can be interpreted as accounting for the emission (absorption) 
by the phonons of m (m') quanta with frequency uj q . Since we are interested in the net number of quanta absorption 
by the phonon bath n = m! — m, it is natural to define 

oo 

C n )P h(r) ^ E - ™> »)<W,phM ( A3 ) 

m,m' 

with S(m, m r ) denoting the Kroncckcr symbol. The identity 

f(a,b,n) = V — — — -5{m! -m,n) 



mim' 

m.m'—0 



^—^ ml (m+n)! — 
m— 

^ r , fc '" + " : n < 

m! (m+n)! 




m— —n 
/2 



J n (2v^h) , (A4) 



where J n (x) denotes the modified Bessel function of the first kind, takes the summation boundaries into account 
properly and implies 



C7ph(r) = E Cn.PhM 

aT _N!(i +2nS ) /i + „! 



n— — oo 
+oo 



E 
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U! 



^^nUl + nl)) . (A5) 



For the full phonon correlation function, we can separate the r-dependence as 

-%£(l+2„*) 



n q=l 



E e " m Xh- (A6) 
n 



where = X) • • • 12 such that each phonon mode has a different summation index and n = (m, . . . , ^q), 

71 rti- — oo tlq — — oo 

f2 = (wi, . . . ,u>q) with n ■ J~J = X)^Li n q UJ q- This enables one to express the Fourier transform of the total (electron 
and phonon) bath correlation function 

= / drCr l (r)C ph (r)e +i ^ 



E Mi / ^ e V)e- iwV 



+i(u—n a ■H)Ts~tTl a 

e v ; ° P h 



= £ 7?' °' - "a ■ «) Ch° = E («) (A7) 

as a weighted superposition of purely electronic Fourier transforms that are evaluated at phonon-shifted frequencies. 
Each term 7" n in the sum accounts for a single electronic jump across junction a triggering the net absorption by 
the phonons of (n Q ,i, . . . , n a _o) quanta carrying a net energy of e a — n a ■ Q. In particular, for a single phonon mode 
(Q = 1) this leads to Eq. (5) in the manuscript. 

The discussion in the manuscript however also applies for multiple phonon modes at the same temperature: The 
bosonic symmetry relation n B = e~P~ phU ' q (1 + rig) for example implies for the weight factors in Eq. (A6) 

c -n a=e -^n a .n c +n^ (A8) 
which can be used to show that the single-mode KMS relation (6) in the manuscript straightforwardly generalizes to 

72> I+ n a (-w) = e-^-^ n ^e^ n ^^ h _ na {+uj), (A9) 

where e a = n a ■ $7 is the total net energy absorbed by the phonon bath when an electronic transition with lead 
a occurs. For the total entropy production, this total transferred energy e Q is relevant. Therefore, when inserting 
counting fields in the Liouvillian it suffices to count the net bosonic energy transfer triggered by electronic jumps 
across each junction instead of counting each phonon mode separately, such that Eq. (7) in the manuscript generalizes 
to 

= E E (aio) 

a£{L,R} n a 

-7?2,n Q (-?) l2i.n a (+e)e +iXo e+ ie ^ 

-7i\n Q e+ ie °«" - 72 \n a (+e) 

Q l h |2 

with the renormalized dot level e = e — Y q . The dimensionless electronic counting fields y Q enable one to extract 

9=1 

the complete electronic particle counting statistics. In contrast, the phonon energy counting fields £ Q have dimension 
of inverse energy and account for the statistics of phonon-related energy transfers triggered by electronic jumps across 
junction a. For the case of a single phonon mode (see main manuscript), the total phonon energy and phonon number 
for electronic jumps across junction a are tightly coupled e a = n a fl, such that the full energy counting statistics 
for the phonons also yields the full particle counting statistics. To obtain the full particle counting statistics for 
multi-mode phonons one would have to make the counting fields mode-dependent e in a^£a — > e 1 H,=i «<*,<,€<*, ^ which 
would exceed the focus of this paper. 
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Appendix B: Full Counting Statistics and the Fluctuation Theorem 

The counting-field dependent Liouvillian (A10) and its single-mode version (7) enable for the construction of the 
FCS, i.e., the probability P (l) (z/> (R)(i) to observe net electrons transfers to lead a and a net associated 

™ol > n el ' e ph ' e ph 

energy absorbed by the phonon reservoir. We are interested in the large-time limit, due to charge conservation, 
it suffices to consider a single electronic counting field. In the following, we will therefore consider xl = and 
Xr = X leading to the probability distribution P (i ) (n){t) associated with nj^ net electrons emitted into the 

right reservoir. It is technically convenient to construct cumulants of the probability distribution via the cumulant- 
generating function (CGF), which in the long-term limit becomes 

C&.&.&.f)'^ (Bl) 

where A(x,£l,£r) denotes the dominant eigenvalue of the Liouvillian which vanishes when all counting fields are set 
to zero A(0, 0, 0) = 0. Cumulants are obtained by taking derivatives with respect to the counting field of interest: The 
trivial long-term time-dependence then enables one to consider the CGF for the currents A(x,£l,£r) instead, such 
that the electronic particle current and noise at the right junction for example may be obtained via 

I = (-i)d x A( x ,0,0)| x=0 , S = (-i) 2 3l A( X ,0,0)| X=0 . (B2) 

Similarly, the total phonon energy current and noise follow by disregarding the difference between left- and right- 
associated emissions 

/ P h = A((U,OI C=0 . ^ = (-i) 2 9| A((U,Ol f=0 ■ ( B3 ) 

The full distribution may be obtained from the CGF by performing an inverse Fourier transform, which in the 
long-term limit becomes 

+ 7T +00 +00 



el i> '""ph 



From properties of the Fourier transform it follows that a shift-symmetry of the CGF is associated with a fluctuation 
theorem for the FCS, i.e., 

A(-X,-£l,-6i) = \(+X + iAHL+iA L ,+£R + iA R ) 

V 1 +"ci '+ e ph '+ e P h (R) A | (L) A . (R) A mrx 

bmln- — = n^'A + e^AL+e^AR. (B5) 

el ' ph ' ph 

We aim at obtaining the affinities A, Al, and Ar, which with Eq. (Bl) yield a long-term symmetry of the CGF. In 
our model (A10), it is most convenient to derive the corresponding symmetry from the characteristic polynomial of 
the Liouvillian. 

To see it we first write the off-diagonal matrix elements of the Liouvillian as explicit functions of the counting fields 
(the diagonal entries do not depend on the counting fields) 

C 12 (+X,Hl,Hr) = E(^i,n(+f)e +i " na +72 i l ) n(+e)e +ix e +i " n ^) , 

n 

C 2 i(+X, Hl, Hr) = £ (7i L 2,„M)e +i '^ + 7f 2 ,„M)e- ix e +i "^) , (B6) 

n 

where we recall that xl = and xr = X- Inserting the KMS condition (A9) and replacing n — > — n under the sum 
we obtain 



■ e 



+7il ( — e)e~ ix e +inn ^ R e~^ pll ~P R ^ nn e l3R ^~^ n ^ 

e+W-^Cn (x + LA, £ L + iA L ,$ R + iA R ) (B7) 
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with A = f3 R (e — (Ar) — j8i(e — /Ui), = /3 p h — Pl, and A# = /3 p h — /3r. This symmetry is equivalent to 

C2i(-X,-^L,-U) = e-^-^)£ 12 ( x + iA,Z L +iA L ,U + iA R ). (B8) 

The characteristic polynomial can be written as 

V(X,tL,tR) = \C( X ^l^r) - Xl\ 

= (Ai-A)(£ 22 -A)-£ 21 ( X ,a,^)A 2 (x,a,^) I (B9) 

where we see that when evaluating the expression at negative arguments the exponential prefactors from Eqns. (B7) 
and (B8) cancel and we easily read off the symmetry 

V(-x, -60 = V(+X + iA, Hl + LA L ,+£ R + iA R ) . (BIO) 

Since the two eigenvalues of the Liouvillian A(x, £lj£r) and A(x, £,l,£,r) are given by the roots of the characteristic 
polynomial T>(x,£,l,^r) = [A — A(x,£l,£r)] [A — A(x, £l, > this symmetry transfers to the long-term CGF as 

K-X,-Zl,-Zr) = \(+X + iAHL+iA L ,+tR + iA R ) (Bll) 

and eventually implies the validity of the full fluctuation theorem (8) in the main text. 



